function v = velocity_centre_tri(p,x,K,darcy_para)
% M = [ones(3,1) x];
% Mi = inv(M);
% dp = Mi(2:3,:)*p;
% 
% v = -1/(darcy_para.mu*darcy_para.phi)*([K(1) K(3);K(3) K(2)]*dp);

% x(1,1)*x(2,2) - x(2,1)*x(1,2) - x(1,1)*x(3,2) + x(3,1)*x(1,2) + x(2,1)*x(3,2) - x(3,1)*x(2,2)
detM = x(1,1)*x(2,2) - x(2,1)*x(1,2) - x(1,1)*x(3,2) + x(3,1)*x(1,2) + x(2,1)*x(3,2) - x(3,1)*x(2,2);
 
v = -1/(detM*darcy_para.mu*darcy_para.phi)*...
    [K(1)*(p(1)*(x(2,2) - x(3,2)) - p(2)*(x(1,2) - x(3,2)) + p(3)*(x(1,2) - x(2,2))) - K(3)*(p(1)*(x(2,1) - x(3,1)) - p(2)*(x(1,1) - x(3,1)) + p(3)*(x(1,1) - x(2,1)));...
     K(3)*(p(1)*(x(2,2) - x(3,2)) - p(2)*(x(1,2) - x(3,2)) + p(3)*(x(1,2) - x(2,2))) - K(2)*(p(1)*(x(2,1) - x(3,1)) - p(2)*(x(1,1) - x(3,1)) + p(3)*(x(1,1) - x(2,1)))];
 
